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Abstract. A two-species spatially extended system of hosts and parasitoids is 
studied. There are two distinct kinds of coexistence; one with populations distributed 
homogeneously in space and another one with spatiotemporal patterns. In the latter 
case, there are noise-sustained oscillations in the population densities, whereas in the 
former one the densities are essentially constants in time with small fluctuations. We 
introduce several metrics to characterize the patterns and onset thereof. We also build 
a consistent sequence of corrections to the mean-field equations using a posteriori 
knowledge from simulations. These corrections both lead to better description of the 
dynamics and connect the patterns to it. The analysis is readily applicable to realistic 
systems, which we demonstrate by an example using an empirical metapopulation 
landscape. 
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1. Introduction 

Pattern formation has gained a lot of interest in both the physical [IIEIEIIIIEIEIITIIEI 
U EE EH El [131 M, QUI US] and the ecological [TZl HSl HSl 1201 EH 1221 literature. Patterns 
may either emerge from the dynamics itself or reflect the structure of the underlying 
landscape. Since the mechanism of pattern formation plays a crucial role, for instance, 
in ecological systems when one is concerned about extinctions of species [TTJ [23j [21], it 
is of fundamental interest to be able to pin down the essential mechanism. 

The classical model of two-species systems, such as predator-prey or host-parasitoid 
systems, is the Lotka-Volterra (LV) model [25]. It assumes fully stirred, or panmictic, 
species. However, introducing explicit space typically leads to correlated structures 
where the assumption breaks down [IT]. The correlations manifest themselves in a 
wide variety of forms. The most common ones are spray-like or flame-like in spatial 
dimensions [H [191 1201 122] , ripple-like in space-time [20], and spiral-like [To] [T8l I2T1 EM] . 
Similar patterning has also been observed in several related models in statistical physics 
[31 [101 El [121 UHL i n chemical surface catalysis [21 SI El [Zl EEE1 E], in calcium signaling 
in cells [26l [271 EE], and in the infamous complex Ginzburg-Landau equation (CGLE) 
[29] . for instance. 

In ecological and chemical systems, the correlations typically weaken the 
interactions, since species tend to be aggregated within themselves. They also 
provide the prey (host) a spatial refuge since around the prey (host) there are less 
predators (parasitoids). This is called the spatial rescue effect [IT]. Altogether, 
spatial inhomogeneity can stabilize the dynamics and strengthen population coexistence 
[XT] [20] [30] even via several different mechanisms [23J . Spatial two-population dynamics 
can also lead to counter-intuitive effects in which increasing the host (prey in prey- 
predator systems) spreading rate actually leads to smaller host population sizes [31 J . 
Understanding spatiotemporal dynamics in highly fragmented landscapes [32] is of even 
greater difficulty 

There are also myriads of empirical observations of the patterning in ecological 
systems. Among the most intriguing examples are field voles in Northern Britain with 
evidence of traveling waves [33], mussel beds in the Wadden Sea in the Netherlands 
with regular spatial patterns [31] together with voles [35] and lemmings [36j in Northern 
Europe. Recently, similar studies have also been performed in experimental laboratory 
conditions [5Tj . 

In this contribution, we study a numerical model of spatially extended two- 
population, or three-state, dynamics. We formulate the model in terms of hosts and 
parasitoids but in the scope of this work, these are interchangeable with prey and 
predator, respectively. It is defined on a regular square lattice, and the spreading of the 
species occurs in a distance- dependent way according to the corresponding incidence 
functions (see below). There are two regimes of coexistence: the species can be either 
distributed homogeneously in space, or build up spatial correlations or patterns (see 
Fig. CD . 
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Figure 1. Two snapshots of the system with coexistence. Left: disordered 
homogeneous structure for parameters Wh = 3.0, w p = 1.5, Xh — 0.63, X p = 1.3, 
and d = 0.9. Right: a patterned state; parameters are as in the non-patterned one 
except for X p — 2.5. Sites with e are white, gray stands for h, and p is shown in black. 

We study the system from two points of view. The first one concentrates on the 
patterns, which we first coarse-grain into areas dominated by one of the species or empty 
space. They give rise to vortices as the corner points of three different areas, and domain 
walls as the boundary lines between them. We introduce several quantities describing 
the geometry and dynamics of these objects. They include the instantaneous velocities 
of the vortices, their number and lifetime, the average length of the domain walls and 
the shape of the dominance areas close to the vortices. All these quantities support the 
division of the parameter space into two kinds of coexistence. 

The second point of view is that of the global dynamics of the system. Here, we 
extend on our previous results [38J in which we have shown that for large enough systems 
there are no non-linearities and in the spatially correlated coexistence regime, the global 
dynamics is governed by noise-sustained oscillations not conforming to the traditional 
limit cycles. See Fig. |2] for examples of the corresponding time series. We show that 
the linearization coefficients differ from what one gets from the complete mixing, or 
mean-field, assumption, that this is caused by the dependence of the effective spreading 
parameters on the instantaneous global population densities, and that this dependence 
is vastly different in the two different coexistence regimes. We also link the dynamics 
together with the patterns by showing that the former one can be understood as an 
infinite stream of aperiodically occurring spontaneous synchronizations. 

Both the tools for the patterns classification and the analysis of the dynamics are 
not restricted to the particular case studied here. Foremost, they are not restricted 
to two species (or three states) or two spatial dimensions. In addition, their domain 
of applicability does not lie in ecology only; the voter models of statistical physics 
PE El dOl El EEl [15] and certain surface catalysis reactions [2, HI El 0, El E] lead to similar 
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Figure 2. Upper panel: the population densities in the patterned state for a system 
of size L — 512 (solid lines), and a subsystem of size L = 64 (dashed lines) of a similar 
one. The upper curves correspond to hosts and the lower ones to parasitoids. The 
parameters are as in Fig.[TJ Lower panel: the population densities in the homogeneous 
state for the hosts (solid line) and parasitoids (dashed line) for a system of size L = 512. 
The latter line has been shifted upwards for clarity. The parameters are as in the upper 
panel except for A p = 1.3. 



patterning, and sometimes similar dynamics, as well, and the machinery discussed here 
provides means to characterize and classify the patterns also in these areas. 

This paper is organized as follows. In Sec. II the model is defined in detail together 
with the coarse-graining procedure, vortices, domain walls and quantities derived from 
them. The section also discusses the program to correct the MF equations. In Sees. Ill 
and IV the results regarding patterns and dynamics, respectively, are presented and 
discussed. Sec. V contains an example application of the analysis, and Sec. VI discusses 
the results both from the physical and the ecological point of view and makes connections 
to earlier and contemporary literature. Finally, Sec. VII wraps up the paper and 
concludes. 

2. Model 

2.1. Definition 

The model describes annual host-par asitoid dynamics on a regular two-dimensional 
square lattice A. It is inspired by Ref. [39] but has a wider interaction range as typical 
in metapopulation dynamics [32]. At a given time, a lattice site can be empty (state e), 
populated by a host either without (h) or with (p) parasitoids. The dynamics allows for 
a cycle of transitions, e^/i— >e— > . . .. This is a simplification that neglects the 
decay of the hosts on their own, and in turn emphasizes the role of the parasitoids. In 
particular, spontaneous deaths of non-parasitized hosts are assumed to be rare enough to 
be considered nonexistent. Even though this means that hosts live forever if parasitoids 
are absent, the restriction is not serious. In coexistence it reduces to assuming faster 
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typical extinction times for the parasitoids than for the hosts. An equivalent description 
of the model is as a variant of the susceptible-infected-recovered (SIR) model [ID] with 
longer spreading lengths and rebirth. Here, the infected state corresponds to hosts, the 
recovered one to parasitoids, and the rebirth to the spontaneous death of the parasitoids. 

The transition probabilities between the states are functions of the connectivities 
which in turn depend on the local populations and are thus directly related to the number 
of immigrants arriving at a lattice site. For a given configuration, the connectivity of 
species h or p on lattice site x at time t is 

4| p (x, t) = J2 ^b(l x - x 'l) t) (1) 

x' 

where x/i|p( x ;^) is the characteristic function, i.e. = 1 if the state of x at time t is h or 
p, respectively, and = else. The kernel has an exponential decay with a characteristic 
scale Wh\ p 

fc A |p(|x-x'|) ocexpf-— — —J, (2) 

and is normalized such that its integral over the whole plane equals one. In time step 
t — > t + 1, the transition e — > h takes place with probability Xhlh and transition h —>■ p 
with probability X P I P . The parasitoids die out (the transition p — > e) with probability 
a irrespective of the surroundings. Parallel updates are used. There is an absorbing 
state with the lattice completely filled with hosts without any parasitoids, together with 
long-lived reactive coexistence states. 

There are numerous assumptions in the model from the ecological point of 
view. The first and foremost are the long-range interactions. In realistic cases, they 
are particularly important for the stability of the system since they contribute to 
immigration and survival also in remote habitat patches. In a regular lattice, their 
role is not as emphasized. The assumption, however, is fairly weak, since the dispersal 
takes place according to an exponentially decaying kernel, i.e. there is a typical short 
interaction length. Furthermore, the connectivity (JTJ) is defined such that parasitized 
hosts do not contribute to the host connectivity. Naturally, this is not completely true 
in empirical systems, since unparasitized hosts can be found also in parasitized habitat 
patches and these are fully capable of contributing to the spreading of the hosts. All 
these assumptions are minor ones and they do not hinder one from being able to tackle 
the relevant and interesting properties of the model system. 

In the MF approximation, the behavior of the system can be calculated as follows. 
The rate equations for the population densities are 

ht+i = h t + k(1 -h t - p t )h t - np t h t 

Pt+i =Pt~ <5p t + fipth , (3) 

where the coupling constants k and \i are just the total transition rates obtained under 
constant densities. There are three fixed points: the coexistence 

h — - and p = — f- (4) 

+ jx) 
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if /i = fi(h,p) > 5, the extinction of parasitoids 

h — 1 and p = (5) 
otherwise, and the extinction of both species 

h = and p = (6) 

which is achieved if the hosts die out before the parasitoids do. The coexistence fixed 
point can be either stable or unstable depending on the parameters. The unstable case 
corresponds to a limit cycle. Below, we compare these elementary observations to the 
behavior of the spatially extended system. Since the MF approximation is merely a 
rough one, there is no particular reason to expect one-to-one correspondence. 



2. 2. Characterizing patterns 

To characterize the patterns, consider a spatially smoothed continuous field of 
population densities [38] 

Ph\ P ,w(^ *) = Ml x - x 'l) Xh\ P (x', t) , (7) 

x' 

where the smoothing kernel k w (x) is a two-dimensional Gaussian such that X^x k w {x) = 
1, X)x x ^™( x ) = 0, and J]x x2 ^( x ) = VJ% - Here the variance w is called the smoothing 
width, and the summation runs over all lattice sites. At each site x, the smoothed 
densities ph\ PtW (x.,t) oscillate around the temporally and spatially averaged densities 

(8) 

and 

1 1 T 

We use these to divide the two-dimensional phase space spanned by the densities into 
three sectors, defined in the caption of Fig. [31 They are then used to divide the system 
into prevalence domains so that site x at time t is defined to belong to a domain 
according to the region of the phase space containing the smoothed densities at location 
x, i.e. (ph(x, t), p p (x, £)) in Eq. (JJJ). The regions coarse-grain on a scale somewhat larger 
than the typical interaction lengths. They are not sensitive to changes in the smoothing 
width w given that it is larger than the interaction lengths and smaller than the system 
size L. 

The domains readily give rise to more sophisticated quantities. To start, define 
the vortices as the corner points of the three different types of domains. These are 
associated with a positive or negative unit "charge" since one can encounter the three 
domains in two different orders by traversing a circle around it in a given direction (say, 
counter-clockwise). Similar structures with three kinds of domains rotating around 
their corner points without any coarse-graining have been identified earlier in statistical 
physics in the context of the three-state voter model [5] , the extended three-state voter 
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Figure 3. (a) The division of the phase space into three sectors associated with 
the three states, e, h and p. For given average densities h and p, the first quadrant 
of the (ph, p p )-plane is divided into three regions by three lines. That separating h- 
and p-dominated regions starts from the average (h,p) (the black dot), goes towards 
increasing ph and p p and is such that its continuation (the dotted line) passes through 
the origin. The other two lines form 120-degree angles with the first one and each 
other, (b) The computed dominance regions of the configuration shown in Fig. [TJ 

model and the three-state Potts model [13] , and combinations thereof [TQl [15] . Also a 
four-state model with game-theoretical inspirations [H] shows similar vortices. In three 
dimensions, the vortices generalize to strings pp. Further, the domain walls are defined 
as the boundary lines between two domains of different types. There are three kinds of 
walls, according to the possible three pairs of domain types they separate, and exactly 
one wall of any given type emanates from a given vortex. The vortices and domain walls 
provide simple means to describe the population patterns, as we demonstrate below. 

To derive quantities that follow the evolution of the vortices and domain walls in 
time, it is necessary to be able to reliably track their motion. To this end, we have 
implemented the following procedure. Let the A_ and A + be the sets of vortices with 
negative and positive sign, respectively, at time t. Denote by B_ and B + the same 
after one discrete time step, i.e. at time t + 1. An interpretation for the movements of 
the vortices, including creation and annihilation or pairs, is a simple pairing of the sets 
X = A_UB + and Y = A + UB_. If an element (x, y) of this pairing consists of vortices in 
A_ and B- (A + and B + ) a negatively (positively) charged vortex is interpreted to have 
moved but not annihilated. If the element consists of vortices in A_ and A + , the two 
vortices are considered to have annihilated with each other, and, finally, if the element 
is a pair of vortices from B_ and B + , the vortices are considered to be newly created 
vortices at time step t + 1 that did not exist at time t. 

There is an enormous number of such pairings, and physically motivated means 
of selecting one are needed. To this end, define a fully connected weighted bipartite 
graph G = (X,Y, E,w), where the vortex sets X and Y are the vertices, the edge set 
their Cartesian product E = X x Y , and the weight function w assigns each edge (x, y) 
a weight equal to the Euclidean distance between vortices x and y. We use the so- 
called Hungarian method from graph theory to obtain the pairing of the sets X and Y 
that has the minimum total weight. This pairing is then used to interpret the motion, 
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annihilation, and creation of the vortices as described above. The Hungarian method 
itself is a well-known exact optimization algorithm for the bipartite pairing problem 
based on graph flows. It is fairly complex but it runs in cubic time with respect to the 
number of elements in the sets to be paired. Detailed explanations of it can be found 
in the literature [ill 132] . 

Building on the tracking procedure for the vortices, it is simple to outline a similar 
procedure for the domain walls. We consider a given wall observed at time t to be the 
same wall as another one at time t + 1 if and only if it is of the same type (separates 
the same two domain types) and if both its endpoint vortices are identified as the same 
according to the method above. A necessary condition for this is that the vortices have 
not been annihilated. 

After implementing these identification procedures, the vortices and walls can be 
used to derive a practically unlimited number of quantities that, in a way or another, 
describe the dynamics or the statics of the patterns. Here, we have used several, selecting 
the ones considered most suitable. These include the jump length of the vortices, 
i.e. distance traveled by a vortex in unit time, the number and lifetime of vortices, 
the geometry of the domains in the immediate vicinity of the vortices described by the 
widths of the sectors the domains form around them, and the lengths of the domain 
walls themselves. See Fig. H] for a schematic illustration of some of these quantities. 




Figure 4. An illustration of the vortices, the domain walls, and the tangents of them 
describing the geometry of the walls near the vortices. The background is a magnified 
portion of Fig. [3] with all three domains colored white and marked with capital letters, 
and the boundaries marked gray. The locations of the vortices are shown with thick 
black circles. The smaller black circles show lattice sites on the wall that are at the 
distance of wt = 5 lattice units from the vortices, and the dashed lines drawn via the 
vortices and these points are the domain wall tangents. 

Upon applying these metrics, attention must be paid to the fact that random 
fluctuations around the average densities also create vortices and domain walls. This 
is an unfortunate side effect of the definitions that must not be neglected but has to 
be taken into account. This is especially important in the case of the domain wall 
lengths, in which the fluctuations give rise to a nonzero background level, and it is the 
deviations of the wall length statistics from this background - not from zero - that 
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signals the onset of patterns. Here, we take the background into account by creating 
random configurations of the system by assigning each lattice site a state of e, h, or 
p randomly using the spatiotemporally averaged densities h and p from a simulation 
as the occurrence probabilities of the respective states. The background vortex and 
wall length levels are then obtained by computing the same metrics for these random 
configurations. 



2.3. Effective mean-field dynamics 

The MF equations ([3]) discussed in the Introduction assume full mixing both within a 
single species and between two species. In spatially extended systems, such as the 
one studied here, the assumptions fail almost invariably. There are several known 
ways to incorporate spatial effects. These include, for example, pair approximations 
[131 HU HS1 HH1 H7] , rescaling the coordination number of the lattice [M], and others [49J. 
Here, our approach is to use input from numerical simulations to introduce corrections 
of different order to the MF equations. These corrections are such that there is no need 
to specify a functional form for the interactions before the analysis. 

The zeroth-order correction is to consider the change in the effective values of the 
spreading rate parameters introduced by the spatial effects. This is done by running 
the simulations for different values of the parameters, obtaining the average population 
densities as a function of them, and subsequently using the inverse of Eqs. (J4]) to arrive 
at those effective values of the MF parameters that lead to the densities observed in the 
simulations. However, common experience suggests that this kind of rescaling of the 
parameters is an utterly implausible explanation for the spatial effects. Instead, higher 
order schemes have to be used. 

The first-order correction builds on the zeroth-order one. Now, in addition, the 
effective parameters are considered functions of the instantaneous population densities 
h t and p t . For A a J a (x, f) small, they equal the spreading probabilities, and thus dynamics 
can be written as [38] 

ht+i =fh +J2[\ h k h (x)C eh (x,t) - \ p k p (x)C hp (x,t) (10) 
xeA 

and 

p t+ i = (l-S)pt + X p M x ) C hp (x,t) , (11) 

xGA 

where the influence of the connectivities on the prevalence dynamics is expressed by the 
correlation functions 

C a/3 (x, t) = J2 Xa(x', t) X/9 (x + x', t) . (12) 

An approximation of these equations can be written as (cf. Eqs. (13~1)) 
h t+1 = h t + n(h t ,Pt) (1 - ht -p t ) h t - n(h t ,Pt)ptht 

Pt+i =Pt~ Sp t + n(h t ,p t )p t h t . (13) 
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This is an approximation of the usual MF form with the interaction parameters n(h,p) 
and fi(h,p) generalized to be arbitrary functions of the instantaneous densities. The full 
form of them seems to be in general not obtainable analytically. In general, they can 
be nonlinear, and in particular they do not have to conform to the standard LV model 
nor to any ad-hoc approximations discussed in the literature, for example turning the 
coordination number of the lattice into two effective coordination numbers that serve 
as phenomenological parameters separately for the hosts and the parasitoids [48] . To 
arrive at the first-order correction, perform a series expansion of Eqs. ( TTTj) around the 
fixed point (h,p) and introduce auxiliary variables rjt and 7r t such that h t = h + rj t and 
Pt = p + TT t - Doing this, we arrive to linear order at 

Vt+i \ = ( ah,h a h>p \ f rjt \ ^ 
T^t+i J \ a P ,h a P , P J V Pit J 
with the matrix elements 

= l+K—2Kh—(K+fi)p+dhK h(l—h—p) — dhp: hp 
a h, P = — (K+[/,)h — d p K h(l — h—p) — d P n hp 

a p ,h = VP + d h fJ* hp (15) 
a PtP = 1 — 5 + jih + dp/jb hp , 

in which /x, k, their derivatives, and the population densities are evaluated at the fixed 
point (h,p). 

If the four derivatives are set to zero in Eq. (fT5]) . Eqs. (fl4"|) and (|T5l) fall back to the 
mean-field approximation, and the matrix in Eq. (TT4"]) is directly the linearization matrix 
usable for the standard stability analysis. On the other hand, if the derivatives are 
nonzero, the mean-field analysis is not valid anymore. This situation can be interpreted 
as instantaneous densities affecting the spreading rates. 

As we have demonstrated earlier [38] and will elaborate on below, the linear form of 
Eqs. ( fl4l) and ( fT5l) can be readily fitted to temporal data from simulations of the spatially 
extended model, and numerical values for the matrix elements can be easily obtained. 
This allows one to solve for the four derivatives, k and \i with respect to h and p, and 
subsequently arrive at the first order correction to the MF equations. Note that this 
computation reveals the domain of applicability of the mean-field treatment by either 
yielding vanishing derivatives (applicable) or non- vanishing ones (not applicable). In 
principle, the procedure could be continued ad infinitum to higher derivatives leading to 
higher-order corrections. These calculations are omitted here, however, since the linear 
treatment turns out to be enough to demonstrate and discuss the effect of the corrections. 
Note that the program outlined here is by no means restricted to the present model, 
but it can be applied to virtually any other one with a similar structure. In particular, 
it is not restricted to two-species (or three-state) models, two spatial dimensions and 
not even to discrete time, since continuous-time results can always be retrospectively 
discretized by sampling snapshots, nor discrete space. Promising candidates for such an 
application include those in Refs. [3, El [10l El H21 [15l [50l EH [52l [53] , for example. 
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Figure 5. Distribution of the jump length of the vortices, i.e. the distance traveled 
by a vortex in unit time, in both the non-patterned (A p = 1.3) and the patterned 
(Xp = 2.5) case. Other parameters are as in Fig. [TJ The average jump lengths are 
approximately 3.6 and 9.2 lattice units for X p — 2.5 and A p = 1.3, respectively. 
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Figure 6. Cumulative distribution of the vortex lifetime according to the tracking 
algorithm in non-patterned and patterned cases. See the caption of Fig. [5] for the 
values of the parameters. The average lifetimes are approximately 2.9 and 8.7 discrete 
time units for A p = 1.3 and A p = 2.5, respectively. 

The vortices and the domain walls can be located and followed through time as 
described in Section 12.21 In general, considering quantities derived from these leads to 
an observation of two distinct kinds of coexistence regimes. There is one with formation 
of chaotic moving fronts, patterns, and another one without them. All computations 
discussed in detail below support the division. 

First, we consider the jump length of the vortices, i.e. the distance a given vortex 
moves in unit time. A histogram of these is plotted in Fig. [5] for the values of the 
parasitoid spreading rate parameter X p corresponding to the patterned and the non- 
patterned cases. Both distributions are exponential, but the patterned case has a clearly 
smaller average jump length than the non-patterned one. Similarly, we plot the lifetime 
distributions of the vortices in the same two cases in Fig. [6j Now, the difference between 
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them is even bigger; in the patterned case the vortices live typically ten times as long as 
in the non-patterned case. Therefore, with patterns the vortices are more stable than 
without them both with respect to movements and annihilation. Were the typical vortex 
jump lengths considered in units of the typical pattern size (for example, the average 
domain wall length discussed below), the difference between the two cases would be 
even more pronounced. 

In a more detailed inspection, the vortex lifetime distribution in the patterned 
case (A p = 2.5) in Fig. [6] is composed of two parts; for short lifetimes the distribution 
follows that of the non-patterned case, and for long lifetimes, there is a broad tail. A 
direct consideration of the definition of a vortex reveals the cause. Namely, vortices can, 
roughly speaking, originate in two different ways. They can be either signs of a "real" 
inhomogeneous structure such as seen in Fig. (TJ or merely random fluctuations around 
the average population densities. The former leads to the wide tail, whereas the latter 
shows up as the short-lifetime part of the distribution. 
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Figure 7. The distributions of the angles between the tangent directions of the three 
domain walls at a vortex (see Fig. UJ. The angles are multiplied by the vortex charge 
so that for each vortex the three domains are passed in the same direction. The data 
is from the patterned state, and parameters are as in Fig. [1] 



Next, we turn to the geometry of the dominance regions in the vicinity of the 
vortices. We compute the angles between the tangent lines of the three walls (see 
Fig. HJ) and plot the distributions in the patterned case in Fig. [7J The parasite-dominated 
sectors appear to be rather narrow (typically 60 degrees) in comparison to the host- and 
empty-dominated ones, whose typical spread equals approximately 160 degrees. This 
is somewhat balanced by the larger variance of the width of the parasite-dominated 
sectors. In any case, these observations conform to those one can make by studying the 
appearance of Fig. [3] in the vicinity of the vortices by bare eye. 

The ratio of the measured domain wall length and its counterpart in random 
configurations is plotted in Fig. [8] as a function of the two spreading rate parameters, A^ 
and X p . The parameter space is clearly divided into two regions with a continuous 
crossover in between by this quantity. In one of these parts, the ratio equals 
approximately one (corresponding to homogeneous population densities), whereas in the 
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Figure 8. The ratio of the measured average domain wall length to its counterpart 
in random configurations as a function of the parasitoid and host spreading rate 
parameters X p and A/j . The ratio is approximately one in the uniformly distributed case 
and increases to approximately three in the patterned case. In all cases, the distribution 
of the domain wall lengths decays exponentially (not shown). Other parameters are 
as in Fig. [TJ 

other, the ratio clearly differs from one (corresponding to pattern formation). In other 
words, in cases where the system does not exhibit patterns visible to human eye, the 
domain wall length statistics also coincides with that in a spatially homogeneous case, 
whereas with patterns there is a significant difference. Note, however, that the numerical 
value of the ratio bears no meaning by itself, since it can be tuned to almost any value 
by altering the smoothing kernel (in particular the smoothing width w) in the coarse- 
graining procedure. Nevertheless, for any usable value of w, the ratio clearly deviates 
from one and therefore is a usable metric of the patterning. The transition between the 
cases is continuous everywhere on the transition line. We have also considered similar 
statistics for the average vortex jump length and the average vortex lifetimes, recovering 
a continuous crossover in each case. 

Given all these results on the vortices, the domain walls, and their dynamics, there 
remains the question of how much of them can be explained by considering the vortices 
as identical random walkers with random signs in a plane. To study this point in detail, 
we have performed numerical experiments consisting of a set of simulated random walks 
for a given value of the simulation parameters as follows. For each walk, draw its 
duration N rw randomly using the measured vortex lifetime distribution. Let each walk 
start at the origin and perform N rw random steps. For each step, draw the jump length 
from the corresponding measured jump length distribution and choose the direction 
in two dimensions randomly from a uniform distribution. Record the average total 
length of the walks for each case of interest, i.e. for different simulation parameters (and 
consequently different measured distributions). 

In the non-patterned and patterned cases, we get average total walk lengths of 20 
and 9 lattice units, respectively, in the random walk experiment. In the non-patterned 
case, this roughly corresponds to half of the average domain wall length which, in 
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turn, approximately equals the average inter-domain distance. Therefore, the picture of 
vortices as merely uncorrelated random walkers cannot be ignored. On the other hand, 
there is more than an order of magnitude of difference between the average domain wall 
and random walk lengths. This rules out any explanation of the dynamics of the system 
that is based on the random walk picture only. 

Finally, we take a look at the behavior of the number of vortices as a function of 
time. A relevant benchmark for comparison is the complex Ginzburg-Landau equation 
(CGLE) which has been studied extensively (for a good review, see [29] ). Both in the 
usual deterministic version of the equation [54] , and in a version with added noise |55j, 
the number of pairs of vortices n(t) has been observed to be a Markov process with 
77,(i)-dependent creation and annihilation rates. An essential feature of such processes is 
the lack of memory of any kind. We have partially repeated the analysis of Refs. [541 [55] 
in the present case. Average creation and annihilation rates as a function of n(t) can be 
measured once a sufficient number of realizations of the process have been simulated. 
Doing this, we arrive at a description similar to Eqs. (3) and (4) of [55] . In other words, 
the creation rate is a constant and the annihilation rate assumes the form 

S_(n) = An 2 + Bn (16) 

where A and B are constants. In this respect the number of vortices behaves here in 
the same way as in the CGLE. However, a detailed look reveals again an important 
difference. In the present case, the number of vortices as a function of time contains 
an oscillatory component clearly visible both in the time series itself and its Fourier 
power spectrum. This cannot be explained by a one-dimensional memoryless dynamical 
system, and therefore the vortex number statistics lies outside the realm of applicability 
of the analysis of Refs. [3U 133] . 

Instead, to recreate the vortex number process, short-term memory has to be 
incorporated in a way or another. One way to do this is to consider a two-dimensional 
process where the role of the second variable is played by the first discrete time derivative 
of the vortex number n. This leads to a description similar to Eqs. (fill) and (fl3|) for 
the vortex number and its derivative. On the other hand, the same process can be 
formulated in terms of the second time derivative or a temporal delay. In any case, 
one-dimensional processes such as Eq. ffT6l) do not suffice. 

4. Dynamics 

4-1- Effective mean-field 

To gain insight on the dynamics of the system, and to produce the effective low-order 
corrections to the MF equations, we have used Poincare maps [56] reconstructed from 
simulations. Here, the dynamical (scalar) variables of the system at time t + 1 are 
considered to be functions of the same at time t. This mapping is then numerically 
computed from measurements of the dynamical variables from the simulations of the 
model. In particular, in the case of two dynamical variables, a graphical representation 
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of the Poincare maps boils down to two three-dimensional scatter plots. We draw 
them for two different system sizes in the patterned case and for comparison in the 
non-patterned one in Fig. [9j Immediate conclusions can be made. There is indeed 
a clear functional dependence of the densities on their previously assumed values. In 
other words, the dynamics can be described by using the densities themselves only; as a 
dynamical system the present one is two-dimensional. This would not have been obvious 
given the 2L 2 discrete degrees of freedom in the spatially extended system. 

We have also verified this conclusion using attractor reconstruction [57j. There, 
a simulated time series is studied using a variable number N d of degrees of freedom, 
which are the time series itself as such, and with suitable delays r, 2r, 3r, and so on. 
Here, we have used a quarter of the time of a period as the time delay r. With two 
dynamical variables, h(t) and h(t — r) we find that there is a unique mapping (up to 
noise) from their values at time t to their values at time t + 1. Such mapping does 
not exist if only one dynamical variable is used, and these findings hold also when p(t) 
(and the corresponding delayed variable) is used instead of h(t). Therefore, attractor 
reconstruction duplicates the conclusion from above that the system has two degrees of 
freedom as a macroscopic dynamical system. 

Based on the numerical observations, the dynamics can be cast as an analogous 
pair of equations to ([3]) but with different functional forms. Even further, the plots are 
clearly linear in both variables for a large system size (here L = 512, Figs. [9^ and [9b), 
whereas for smaller system sizes (L = 64, Figs. [9fc and(9H) there are deviations from the 
linear form. The point clouds are stretched towards the absorbing state of the system 
with extinct parasitoids. The necessary system size for the onset of linearity can be 
crudely estimated by visual inspection. It reveals the rough approximation L\ in ^ 120, 
slightly larger than the corresponding average domain wall length (see Sec. |3]). Thus, 
the linearity is obeyed, if the system is large enough to accommodate several domains. 

To arrive at an effective MF iteration based on the Poincare plots valid for large 
systems, we perform least-squares fits to them and arrive at Eq. (fl4l) with numerically 
determined coefficients a a y. This picture takes space implicitly into account since the 
interaction parameters (or the matrix elements) are automatically suitably adapted 
to the form and parametrization of the manifolds. In addition to this deterministic 
description, the system exhibits noise whose magnitude - both perpendicular and 
tangential to the manifolds - is proportional to the square root of the system size. 
Also, after the pure noise component and the deterministic picture above have been 
removed, one is left with a residue that comes from the approximation made in the 
linearization. However, for large system sizes this residue is small when compared to 
the linear deterministic component and the noise. The numerically fitted coefficients 
a a y come into play below during the discussion of the effective parameters and sustained 
oscillations. 

The predictions of the MF theory of Eqs. ([3]) for the host and parasitoids densities 
compared with those obtained from simulations are shown in Fig. [101 The host density 
ph does not depend on the host spreading rate parameter X h according to the MF theory 
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Figure 9. The Poincare maps illustrating the dynamics of the system, (a), (b): ht+i 
and pt+i as a function of ht and pt for a system in the patterned state (parameters 
as in Fig. la) with size L x L = 512 x 512. (c), (d): The same for a system of size 
64 x 64. (e), (f): As in (a) and (b) but for a system in the homogeneous state, i.e. with 
X p = 1.3. The simulations have been run for 3000 time steps and data was gathered 
for the last 1000. 



and as is seen from the figure, this is almost the case also in simulations. On the other 
hand, the parasitoid density p p depends strongly on both spreading rate parameters Xh 
and X p both in the MF theory and in the simulations. Altogether, Fig. [10] shows that 
the MF approximation works surprisingly well in predicting the species densities. Note 
however, that these comparisons are only for the average densities and do not constitute 
evidence of the validity of the MF theory in other respects, such as the stability of the 
fixed point or the nature of oscillations, if any. 

In order to arrive at equations that can also predict these more complicated issues, 
we proceed in two phases. The first one - here called the zeroth-order correction - is 
to ask which values of the parameters have to be inserted into the MF equations so 
that they give the same population densities as the simulations. In other words, the 
question is what are the effective MF spreading rate parameters as a function of the 
simulational ones. To this end, we have run simulations for a wide range of parameters 
and computed the effective host and parasitoid spreading rates K e ff and /x e fr by using 
the inverse of Eq. The results are plotted in Fig. [TTJ The conclusions are evident. 
The dependence is strongly nonlinear, as was expected, and the immediate observation 
is therefore that the effective parameters cannot be obtained from the simulational ones 
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Figure 10. Prevalences as functions of the spreading rate parameters A^| p , together 
with the MF predictions. The prediction for ph does not depend on A^, whereas the 
for p p predictions for A/ t = 0.625,1.25,1.875,2.5 and 3.125 from bottom to top are 
shown. The color code is the same in both panels. 



by any simple rescaling. Furthermore, the division between the patterned and non- 
patterned regions of the parameter space is clearly visible in Fig. [TTJ It shows up as a 
ridge in K eS and as a bend in fj, e g . 




Figure 11. Effective values of the parameters k (left) and /i (right) as a function of 
the input parameters A^ and X p . 

In addition to what has been discussed above, the effective parameters do depend 
on the population densities within single realizations with fixed parameters. Here, we 
call this dependence the first-order correction, and it is the one that causes changes in 
the effective iteration equations with respect to Eqs. Q. Indeed, given the described 
dependence, k and \i in Eqs. ([3]) are functions of h t and p t as opposed to constants. 
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To compute the first-order correction, linearize Eqs. ([3]) around the reactive fixed point 
assuming that k and fi are functions of the population densities. One arrives at Eq. ( TT4|) . 
in which there are four matrix elements that are linear in the four unknown first 
derivatives of n and ji with respect to h t and p t . The matrix elements are then set 
to equal those obtained from linear fits to the Poincare plots, i.e. the a a y in Eq. ( fT4l) . 
and the unknown derivatives are solved for. 

The results are shown in Fig. [121 i n which for each considered value of the 
parameters Xh and A p the derivatives are drawn as an arrow depicting the magnitude and 
direction of the gradient of k and /i with respect to h t and p t . The auxiliary coordinate 
system on the bottom left corner fixes the orientation of the gradients. The behavior of 
HeS is especially interesting. In the non-patterned case the derivatives practically vanish. 
This means that the rate parameters k and fi are effectively constants and that the MF 
theory is a rather good approximation. This, in turn, means that the populations are 
well-mixed. On the other hand, the derivative of fj, e g with respect to p t clearly deviates 
from zero and assumes negative values in the patterned case. This is a direct sign 
of the insufficiency of the MF theory in this parameter regime. The correlations that 
for a larger number of parasitoids, makes them aggregated within themselves and thus 
decreases the "free boundary" available for spreading constitute one contributing factor 
to this insufficiency 

We have also checked the consistency of the computations above by solving for the 
effective k and fi from the MF iteration equations ([3]) separately for each time step, 
and considered these as a function of the population densities. The results are in an 
excellent numerical agreement with those above. In principle, the fits of the functional 
form in Eq. (JT4l) could be extended to any order, and used together with the higher- 
order counterparts of the analysis of the derivatives of k and fi above to extract the 
derivatives up to any order. 

4-2. Oscillations 

With patterns the population densities oscillate with an amplitude that fluctuates in 
time (Fig. [TBI . Without patterns, the densities fluctuate randomly around a stable 
fixed point. The angular frequency of the oscillation can be measured by performing 
three-dimensional linear fits to the Poincare maps (Fig. [9]) and extracting the imaginary 
part of the eigenvalues of the resulting matrix a a ^i in Eq. ( fl^l) . The MF prediction for 
the frequency is computed similarly. These are compared in Fig. [TH The frequencies 
from the MF theory and the simulations have roughly the same numerical values, but 
the agreement is far from perfect. 

Instead of the oscillations themselves, the most characteristic property of the 
population densities is the fluctuation of the amplitude. The typical amplitude depends 
both on the parameters and the system size (see Fig. [2]). Decreasing the system size 
increases it but keeps the qualitative behavior and the population densities intact. We 
have previously given an explanation for the amplitude fluctuations [38] which we here 
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Figure 12. The first-order dependence of the effective parameters on the population 
densities ph and p p . At each point, the arrow depicts the derivative of K e ff (above) and 
/U e ff (below) according to the auxiliary coordinate system drawn on bottom left. The 
same points of the parameter space are sampled as in Fig. [TT] 
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Figure 13. Time series of host and parasite prevalences, and the synchronization 
order parameter r (Eq. (fT8| ). The parameters are as in Fig.Q] The amplitudes of the 
oscillations correlate strongly with each other and with the values of r. See Fig. [2] for 
the dependence of the amplitudes on the system size. 



briefly rephrase. Namely, both with and without patterns, the fixed point of the effective 
(deterministic) iteration given by the fits to the Poincare maps is stable. As such, there 
should be no oscillations in either case. However, the eigenvalues of the iteration matrix 
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Figure 14. The angular frequency of the sustained oscillations as a function of the 
host spreading rate parameter for several values of A p . The frequency is measured 
as the polar angle of the eigenvalues of the linear fit of the simulation data to Eq. (I14| . 
The measurements are done for a system of size L x L = 512 x 512 which is large 
enough for the linear fits to work well. Other parameters are as in Fig. [TJ The solid, 
dashed, and dotted lines show the frequency predicted by MF theory, i.e. linearization 
of Eq. ^ and numerical solution of the resulting eigensystem. 



are a complex-conjugated pair in both cases, and the transients towards the fixed point 
are oscillatory. Such transients give naturally rise to two independent time scales, one for 
the oscillation (time of a period) and another one for the decay. Denoting the eigenvalue 
pair as pe'^, the ratio of the time scales is 



r osc logp' 

Inspecting this numerically reveals small values without patterns and large values with 
them. Therefore, in the latter case any fluctuations away from the fixed point lead to 
slow oscillatory convergence towards the fixed point which is repeatedly disrupted by 
stochasticity giving rise to a perpetual decay - noise-sustained oscillations. 

Note that when approaching the stability limit from below, p approaches one, and 
the ratio v diverges towards positive infinity. With p > 1, the ratio reassumes finite 
values, now with a negative sign, but does not anymore carry the same physical meaning 
as below the limit. 

The time scale ratio v both measured from simulations and predicted by the MF 
equations and the excess domain wall length (the difference between the average domain 
wall length and its value in random configurations) are plotted in Fig. [15] as a function 
of the parasitoid spreading rate parameter X p . The behavior of v is drastically different 
between simulations and theory. In the latter case, the ratio diverges at around \ = 1.9 
and assumes negative values at higher X p . In other words, the system undergoes a 
bifurcation from a stable fixed point to an unstable one characterized by a limit cycle. 
On the other hand, the simulated time scale ratio v peaks around almost the same point, 
but does not in fact diverge, signaling that there is no limit cycle involved. To see the 
connection between the temporal and spatial properties of the system, note that the 
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dependence of the average domain wall length and the time scale ratio on the parasitoid 
spreading rate parameter X p is similar for A p < 2.2. This can be interpreted as direct 
evidence that the pattern formation and the deviations of the dynamics from the MF 
description are closely intertwined. Also the behavior of the time scale ratio and the 
excess line length is qualitatively similar for X p > 2.2 since both are constants in this 
regime. This is comparable to the coupled discrete oscillators of Wood and co-workers 
[52] [53] . where a related system has a second-order synchronization transition in the 
MF theory and in higher dimensions but not in two. 
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Figure 15. The timescale ratio of Eq. (j 1 T[) in both the simulations and the MF 
theory plotted together with the difference between the average domain wall length in 
the actual simulation and in random configurations. All other parameters except for 
A p are as in Fig. [T] 



There is also another straightforward spatial interpretation of the amplitude 
fluctuations. Consider the system divided into small boxes whose width is clearly smaller 
than the typical width of the coarse-grained patterns, but which at the same time are 
large enough so that local self-averaging occurs inside each box. Now a local phase 
4>i for each box i can be defined as the phase angle of the population densities inside 
the box with respect to the average densities (see Fig. [3]). Given these, the degree of 
synchronization between the boxes can be defined as is customary in the context of the 
Kuramoto model for coupled oscillators [SS] 

r=i|E^I, (18) 
j 

where the summation runs over all N boxes. 

In this computation, we have used the box size I x I = 64 x 64 in accordance with the 
criteria above, and computed the synchronization order parameter r for each time step. 
This is plotted in Fig. [13] together with the host and parasitoid densities in the same 
simulation run. The amplitude of the density oscillations correlates strongly with r. 
This is to be interpreted such that the system spontaneously synchronizes over at least 
intermediate length scales at random intervals. As a feature of secondary importance, 
the synchronization order parameter r oscillates with a time of a period half of that 



Characterizing spatiotemporal patterns in three-state lattice models 



22 



of the population densities, which is also seen in the Fourier power spectrum of r (not 
shown) . 

4-3. Phase diagram 

Given the fact that the noise- sustained oscillations can be interpreted in terms of 
spontaneous synchronization, it is an interesting question whether it is macroscopic 
in character or "merely" a finite- size effect. To address this issue, we have computed 
the temporal average of the synchronization order parameter r (fl8l) for different system 
sizes L both in the patterned and non-patterned cases. The results are shown in Fig. [16j 
In all cases, the scaling r ~ L~ l is found (similarly to [52l [53] ). Thus, we conclude that 
the system does not show long-range order and thus does not completely synchronize 
at macroscopic scales. 
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Figure 16. The synchronization order parameter r of Eq. (fT8|) as a function of system 
size L for several values of X p . In each case, the scaling r ~ L _1 is found, demonstrating 
the lack of long-range order in the system. 

In spite of this, the phase diagram remains interesting. In addition to the extinction 
transition of the parasitoids, the coexistence phase can be split into three. There are 
the non-patterned and the patterned regions, and since the crossover between them is 
not sharp a "gray area" between them. We have built the phase diagram in the (A^, X p )- 
plane by defining the boundary between the non-patterned case and the gray area to be 
the line at which the timescale ratio v of Eq. (TTTj) equals 1.0, and that of the gray area 
and the patterned phase where the ratio equals 4.0. Albeit arbitrary, these definitions 
seem fit since upon looking at the population densities as a function of time and figures 
such as Fig. [H the patterns and the oscillations are clearly not there if the ratio is below 
one, and, on the other hand, if it exceeds 4, both are clearly observed. Note that the 
timescale ratio assumes values as high as 100 (see Fig. [15]). The qualitative features of 
the phase diagrams are not sensitive to altering the limiting values of v. 

The phase diagrams are shown in Fig. [T7| for two choices of the spreading widths 
Wh and w p . It is seen that for any value of Xh the four phases (parasitoid extinction, 
homogeneous coexistence, the gray area, and patterned coexistence) follow each other 
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in this order if the parasite spreading rate parameter X p is increased. This feature does 
not depend on the particular values of the spreading widths, or even their mutual order. 
The phase structure bears some resemblance to related earlier models [IH HHJ H7] in that 
oscillatory and non-oscillatory regions and a boundary between them is recovered. See 
Sec. M for a detailed discussion and comparison. Fig. [TTb also shows the vortex density, 
i.e. number of vortices per lattice site, as a function of X p . One sees that the number of 
vortices decreases rapidly while going from the nonoscillatory region via the gray area 
to the oscillatory regime. 




Figure 17. Phase diagrams with varying Xh and A p , for different spreading ranges, 
(a) Wh = 3 and w p = 1.5, (b) switched Wh = 1-5 and w p — 3. In both cases four 
"phases" are found, parasite extinction, coexistence without oscillations, coexistence 
with oscillations, and a gray area between the two, which is due to the fact the there 
is not, in fact, a sharp transition between the oscillatory and non-oscillatory regimes, 
but a continuous cross-over. The upper and lower bounds of the gray area are defined 
as the lines where the timescale ratio of Eq. ([T"7|) equals 1 and 4, respectively. The 
system size is L x L = 512 x 512. (c) The vortex density as a function of the parasitoid 
spreading rate parameter A p for the case of panel (a) with A/i = 0.63 and the smoothing 
width (7 = 8. 

Finally, we note that even if the system does not in general synchronize 
macroscopically, such synchronization can be artificially achieved by making the system 
small with respect to the spreading widths. This can originate in either making the 
system really small, or keeping the size fixed and increasing the spreading widths. In 
these cases, the MF theory probably works better than in what has been studied here. 
This, however, lies outside the scope of interest since extremely small system sizes or 
very large spreading lengths are needed. 

5. A realistic example 

To give a concrete example of how to apply the methods introduced here, we illustrate 
the role of coherent dynamics in biological systems by resorting to a well-known 
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metapopulation system, investigated at length by the Metapopulation Research Group 
at University of Helsinki, Finland [321 [591 EDI EU E2]- The species in question are the 
Glanville fritillary butterfly Melitaea cinxia and its specialist parasitoid wasp Cotesia 
melitaearum on the network of habitats on the Aland islands in the Baltic Sea (60° 
N, 20° E) between Finland and Sweden. Here we use the known patch geography as 
the "lattice" . The two populations follow a variant of the model used above, with the 
differences that the contribution of each patch to the connectivities is multiplied by 
its area, and that in addition to the connectivity-driven spreading a small amount of 
uniformly random reproduction is used. These differences are biologically motivated. 
Intuitively, the immigration pressure from a given patch increases with the local 
population, which itself is positively coupled to the amount of locally available suitable 
habitat, the patch size. The random spreading mimics occasional long-range dispersal, 
and its effect is to prevent extinction on remote subnetworks, which are naturally present 
at the outskirts of the sparse archipelago. Also, open boundaries are used. This does 
not restrict the validity of the analysis for two reasons. First, the equations (114j) are 
an observation made from the simulations, and by plotting the Poincare maps as in 
Fig. El one sees that the equations hold both for open and periodic boundaries. Second, 
the equations (TT3T) are obtained from Eqs. (fl3l) by straightforward differentiation, and 
Eqs. ( ITT]) , in turn, are a rough approximation equally valid in either case. 

Fig. [18] shows a snapshot of the system, and the inset contains the respective time 
series; see the caption for the parameters. The question is now, do we find evidence of 
oscillations and/or patterns here? A similar analysis of the dynamics, from Eq. ( fl4l) 
is readily carried out, revealing the eigenvalues A± ~ 0.84 ± 0.17i for the parameters 
used in Fig. [TSJ These imply, in accordance with the appearance of the time series, that 
the system contains patterns and exhibits oscillations coming from the combination 
of an unstable fixed point with oscillatory transients and demographic stochasticity. 
The patterning is clearly visible in Fig. [TBI It can be analyzed using the smoothed 
densities by considering the sum in the expression for them, Eq. ([7]), as a summation 
over the patches, and using the position-dependent smoothed densities assuming average 
population densities at each patch as the triple point in the phase space (see Fig. [3^). 
The rest of the analysis can then be carried out as explained in Section I2.2L Doing this, 
one finds a domain structure analogous to that on the regular lattice with the average 
domain wall length (£) /(^random) ~ 1-8. 

These observations are of interest for a couple of reasons. First, the Aland landscape 
is rather heterogeneous, and at a closer inspection (63] appears to consist of a variety 
of densely connected subsystems, which are weakly coupled to each other. Due to the 
differences in the landscape with respect to a regular lattice, it is not clear, a priori, that 
the whole system has similar pattern-included noise-sustained oscillations as the version 
defined on a regular lattice. Instead, the expectation for a single densely connected 
subsystem is to be closer to the mean-field limit [45J obeyed by a fully connected 
(essentially non-spatial) system, and since the number of such dense subsystems is 
relatively small, the expectation is the same for the full system. However, noise- 
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sustained oscillations are observed, and the Aland system works as an example of the 
applicability of the analysis presented here to more realistic cases, than regular lattices. 
This is emphasized by the fact that also the dynamics of the model contains differences 
with respect to the lattice one, in addition to the spatial structure. In any case, the 
eigenvalues of the effective iterative map (Eq. (1141) ) provides information about the 
stability properties and the level of patterning in the system. 
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Figure 18. Main figure: a snapshot of the simulated dynamics of hosts and parasitoids 
on the Aland patch geometry. Each dot corresponds to a single patch. Inset: the 
prevalence of the hosts (dashed line) and parasitoids (solid line) as a function of time. 
As spreading widths we have used Wh — 1000 m, w p — 500 m. The annual death 
probability is 8 = 0.9 and the host and parasitoid spreading rate parameters are 
A/j = 100 and A p = 1400. A fraction of 0.01 of the hosts spreads randomly to all 
patches. The archipelago is roughly 60 km x 80 km in size. 



6. Discussion 

The change between the patterned and non-patterned cases is rather abrupt in terms 
of several quantities. This leads naturally to the question whether the phenomenon is 
actually a phase transition. We argue that this is not the case due to several reasons. 
First, since phase transitions are defined only for infinite systems, the transition in plots 
such as Fig. [8] should become sharper as the system size increases. However the behavior 
remains essentially the same as in Fig. [8] for a wide range of system sizes, except for 
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the fact that for the smallest systems, the highest value the average domain wall length 
assumes becomes smaller due to finite-size effects. In any case, no sharpening of the 
transition is observed, and in this respect the correct description of the behavior is as a 
continuous crossover. 

Second, similar models have been studied with the transition aspect in mind. The 
most relevant one for the present discussion is that by Wood and coworkers |52j [53] , 
which deals with a symmetric three-state model of discrete coupled oscillators. It has 
been found that the system does not have a phase transition in two dimensions but 
does in three and higher dimensions. However, in two dimensions, there is a continuous 
change in the order parameter around the region of the parameter space where the 
transition would take place in higher dimensions. This is similar to what has been 
observed here for the average domain wall lengths, giving further credibility for the 
interpretation of the change from non-oscillatory to oscillatory behavior as a continuous 
cross-over. 

The eigenvalue structure in the non-patterned case also speaks in favor of the 
absence of a sharp transition. We get a complex conjugate eigenvalue pair for the 
matrix a CTj0 -' in Eq. flHj) even in the non-patterned case. The difference lies in the relative 
magnitudes of the oscillation and decay time scales. So, even if the ratio of the timescales 
were used instead of the average domain wall length, there would be no sharp transition 
(cf. Fig. [T5l) . This can be viewed as a consequence of the corrected MF approach 
yielding a stable fixed point both in the patterned and the non-patterned regimes. As 
a conclusion, the difference between the two is strictly speaking only quantitative, or - 
put in other words - in the non-patterned case the spatial correlations are much more 
suppressed than in the patterned one. However, on the other hand, the quantitative 
difference between the cases is rather huge: There are practically no patterns in the non- 
oscillatory systems as can be seen by looking at the system itself (Fig. [1]), the average 
domain wall length (Fig. [8]), and the random walk experiments outlined in Section l2~3l 
for example. Therefore, discussing these two regimes as different appears justified. 

In addition to the discrete-space discrete-time setting studied here, the methods 
of analysis and characterization of patterns and dynamics are equally applicable to 
systems with continuous time or space, for example that in Ref. [49J. Applying the 
methodology to such systems might reveal differences between them and their discrete- 
space counterparts. Studies of these would be interesting already as such. One possible 
cause of the differences is that the discretization of space into lattice cells often comes 
with implicit locally density-dependent establishment, i.e. the local population density is 
restricted such that only one individual can be present at one site at a time. Continuous- 
space approaches typically do not have such restrictions unless explicitly included. 

A very similar system restricted to nearest-neighbour spreading has been recently 
studied [HJ US, H7] both analytically using the pair approximation and numerically. 
The authors of these articles have found that the co-existence phase is divided into 
two regions. These are oscillatory and non-oscillatory coexistence. The similarity to the 
present case is striking. It has been concluded in from the pair approximation 
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that the oscillations are a limit cycle, and indirectly that, since there are oscillations also 
in the spatially extended counterpart, they have to be limit cycles as well. Comparing 
to the present results reveals that such conclusions cannot be made without further 
evidence; the possibility of noise-sustained oscillations remains, and for small systems 
even that of two kinds of oscillatory regimes. While further studies would be needed 
to resolve this issue, note that in some cases the amplitude fluctuates in a way that 
appears similar to the fluctuations in the present case. One example is given by Fig. 9a 
of Ref. On the other hand, there are established results on two-dimensional three- 
state dynamic lattice models with limit cycle oscillations [61]. Furthermore, processes 
on excitable media have been recently characterized in terms of prey-predator systems 
[65J, and the machinery introduced here could provide tools for studying those as well. 

Other comparable approaches taken recently include Refs. j66][67|. In these works, 
a simplified model of two connected patches is studied, and the oscillations are made 
dependent either on the phase or the amplitude in the phase space. Both lead to 
stabilized oscillations. However, Eq. iHM explicitly forbids such dependencies in the 
present case, and therefore the mechanisms suggested in [66l [67] cannot be the root 
cause of the stability in the spatially extended case studied here. 

Another interesting approach is to map related population models [2U IS] to the 
complex Ginzburg-Landau equation |29j. The mappings rely on the existence of an 
unstable fixed point [68] . and the oscillations are classical limit cycles. Therefore, these 
studies and the present one can be described as complementary approaches to each 
other. 

While the present model is without doubt castable as a partial differential equation 
at least in an approximative manner, it is an open question of what the outcome from 
such equations would be. Since the fully-mixed approximation can have both a stable 
and an unstable fixed point, a mapping to the CGLE along the lines of [16] would not 
work, and further complications are also caused by the highly asymmetric reaction rates 
[68] . Due to these reasons, it is not immediately clear that the picture of noise-sustained 
oscillations linked to patterns via spontaneous synchronizations would be captured by 
such a treatment. In any case, this kind of analysis offers an intriguing issue for future 
studies. 

7. Conclusion 

A model of two-species dynamics defined and treated in the language of hosts and 
parasitoids but equally applicable to prey-predator systems as well has been studied. 
In addition to extinct states, the system shows two kinds of coexistence, depending on 
the parameters. These are a patterned state with oscillating population densities and a 
non-patterned state with populations homogeneously distributed in space, and densities 
that fluctuate around their respective averages. This contribution contains two closely 
tied lines of work. One characterizes the patterns via coarse-graining the habitat space 
to domains in which a given species or empty space is dominant. The other concentrates 
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on the dynamics of the system. 

The coarse-grained domains allow for a definition of domain walls as the lines 
separating the domains, and vortices as their corner points. From these, we compute 
several quantities, such as the instantaneous vortex velocities, their lifetime, the lengths 
of the domain walls, and ones describing the geometry of the domains in the vicinity of 
the vortices. Invariably, these quantities reveal the division of the parameter space into 
two kinds of coexistence. We have argued that these two regions are separated by a 
continuous cross-over rather than a conventional phase transition. We have also reasoned 
why the vortex number process is not comparable to that in the CGLE [29l I5H 155] . The 
main argument is based on the fact that the process is not memoryless. 

The dynamics has been studied using Poincare maps [56]. We have found out that as 
a stochastic dynamical system, the present one is two-dimensional. In other words, the 
dynamics is adequately described in terms of the two population densities and the unit- 
time advancing mapping, the Poincare map, of the two. No further dynamical variables 
nor history-dependence are needed. In addition, for large systems the Poincare maps 
are linear leading to an easy characterization of the behavior using eigenvalue analysis. 
In both the patterned and the non-patterned case, this reveals a stable fixed point with 
oscillatory transients. However, between the cases there is a huge difference in the 
level of separation between the associated oscillation and decay time scales. This also 
leads to a direct connection between the patterning and the dynamics, since the ratio 
of the two timescales has been shown to be closely related to the average domain wall 
lengths. The amplitude fluctuations have also been given a spatial interpretation as 
an irregular sequence of spontaneous synchronization and desynchronization events of 
coupled oscillators [58] whose role in this setting is played by mesoscopic subsystems. 

We have also compared the average densities in the system to those predicted by 
a mean-field theory. We have introduced corrections to the MF equations on two levels 
based on series expansions and simulated data. In both zeroth- and first-order settings, 
the cross-over between the patterned and the non-patterned cases is visible and crucial. 

The results have been compared to earlier literature. There are numerous results 
on related systems where oscillations have been recovered. In some of these, the 
oscillations might in fact be noise-sustained in character. Interesting candidates to 
study this possibility include [SJ HEl EH EDI EI]. In all of these, the machinery for 
pattern characterization introduced here could be applied to produce further insight. 

We have applied the analysis briefly to the empirical metapopulation landscape of 
a butterfly and its parasitoids wasp [59], [60] . Future prospects include studying these 
issues further - we have merely highlighted the applicability of the results presented 
here. Studying such systems could, for instance, let one approach the question of how 
to tell the difference between patterns caused by inhomogeneities in the landscape and 
those that are spontaneously formed even on homogeneous substrates. Such extensions 
are naturally motivated also from the ecologists' point of view since metapopulation 
ecology often involves such landscapes [52"] . 
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